home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / toolbox / nmsmax.r < prev    next >
Text File  |  1994-04-25  |  7KB  |  231 lines
  1. #function [x, fmax, nf] = nmsmax(fun, x, stopit, savit)
  2. #NMSMAX  [x, fmax, nf] = NMSMAX(fun, x0, STOPIT, SAVIT) attempts to
  3. #        maximize the function specified by the string fun, using the
  4. #        starting vector x0.  The Nelder-Mead direct search method is used.
  5. #        Output arguments:
  6. #               x    = vector yielding largest function value found,
  7. #               fmax = function value at x,
  8. #               nf   = number of function evaluations.
  9. #        The iteration is terminated when either
  10. #               - the relative size of the simplex is <= STOPIT(1)
  11. #                 (default 1e-3),
  12. #               - STOPIT(2) function evaluations have been performed
  13. #                 (default inf, i.e., no limit), or
  14. #               - a function value equals or exceeds STOPIT(3)
  15. #                 (default inf, i.e., no test on function values).
  16. #        The form of the initial simplex is determined by STOPIT(4):
  17. #           STOPIT(4) = 0: regular simplex (sides of equal length, the default)
  18. #           STOPIT(4) = 1: right-angled simplex.
  19. #        Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
  20. #        If a non-empty fourth parameter string SAVIT is present, then
  21. #        `SAVE SAVIT x fmax nf' is executed after each inner iteration.
  22. #        NB: x0 can be a matrix.  In the output argument, in SAVIT saves,
  23. #            and in function calls, x has the same shape as x0.
  24.  
  25. # References:
  26. # J.E. Dennis, Jr., and D.J. Woods, Optimization on microcomputers:
  27. #     The Nelder-Mead simplex algorithm, in New Computing Environments:
  28. #     Microcomputers in Large-Scale Computing, A. Wouk, ed., Society for
  29. #     Industrial and Applied Mathematics, Philadelphia, 1987, pp. 116-122.
  30. #  N.J. Higham, Optimization by direct search in matrix computations,
  31. #     Numerical Analysis Report No. 197, University of Manchester, UK, 1991;
  32. #     to appear in SIAM J. Matrix Anal. Appl, 14 (2), April 1993.
  33.  
  34. # This is a heavily modified version of FMINS.M supplied with 386-MATLAB
  35. # version 3.5j.
  36.  
  37. # By Nick Higham, Department of Mathematics, University of Manchester, UK.
  38. #                 na.nhigham@na-net.ornl.gov
  39. # July 27, 1991.
  40.  
  41. # Translated to RLaB, Ian Searle
  42. # Feburary 1994.
  43.  
  44. nmsmax = function (fun, X, stopit, savit)
  45. {
  46.   local (beta, vc, ve, fc, scale, fe, x, vk, fmax_old, ...
  47.          fmax, V, trace, fk, vbar, vr, how, tol, size_simplex, ...
  48.          vt, fr, ft, alpha, v1, nf, f, msg, x0, gamma, ...
  49.          j, k, temp, m, n)
  50.  
  51.   x = X;    # Copy the input
  52.   n = prod(size(x));
  53.   x0 = x[:];    # Work with column vector internally.
  54.  
  55.   # Set up convergence parameters etc.
  56.   if (!exist (stopit)) { stopit[1] = 1e-3; }
  57.   tol = stopit[1];    # Tolerance for cgce test based on relative size of simplex.
  58.   if (max(size(stopit)) == 1) { stopit[2] = inf(); }    # Max no. of f-evaluations.
  59.   if (max(size(stopit)) == 2) { stopit[3] = inf(); }    # Default target for f-values.
  60.   if (max(size(stopit)) == 3) { stopit[4] = 0; }    # Default initial simplex.
  61.   if (max(size(stopit)) == 4) { stopit[5] = 1; }    # Default: show progress.
  62.   trace  = stopit[5];
  63.   if (!exist (savit)) { savit = []; }            # File name for snapshots.
  64.  
  65.   V = [zeros(n,1), eye(n,n)];
  66.   f = zeros(n+1,1);
  67.   V[;1] = x0; 
  68.   x = reshape (x0, x.nr, x.nc); 
  69.   f[1] = fun (x);
  70.   fmax_old = f[1];
  71.  
  72.   if (trace) { printf("f(x0) = %9.4e\n", f[1]); }
  73.  
  74.   k = 0; m = 0;
  75.  
  76.   # Set up initial simplex.
  77.   scale = max([norm(x0,"i"),1]);
  78.   if (stopit[4] == 0)
  79.   {
  80.     # Regular simplex - all edges have same length.
  81.     # Generated from construction given in reference [18, pp. 80-81] of [1].
  82.     alpha = scale / (n*sqrt(2)) * [ sqrt(n+1)-1+n, sqrt(n+1)-1 ];
  83.     V[;2:n+1] = (x0 + alpha[2]*ones(n,1)) * ones(1,n);
  84.     for (j in 2:n+1)
  85.     {
  86.       V[j-1;j] = x0[j-1] + alpha[1];
  87.       x = reshape (V[;j], x.nr, x.nc); 
  88.       f[j] = fun (x);
  89.     }
  90.   else
  91.     # Right-angled simplex based on co-ordinate axes.
  92.     alpha = scale*ones(n+1,1);
  93.     for (j in 2:n+1)
  94.     {
  95.       V[;j] = x0 + alpha[j]*V[;j];
  96.       x = reshape (V[;j], x.nr, x.nc); 
  97.       f[j] = fun (x);
  98.     }
  99.   }
  100.   nf = n+1;
  101.   how = "initial  ";
  102.  
  103.   temp = sort (f).val;
  104.   j = sort (f).ind;
  105.   j = j[n+1:1:-1];
  106.   f = f[j]; 
  107.   V = V[;j];
  108.  
  109.   alpha = 1;  beta = 1/2;  gamma = 2;
  110.  
  111.   while (1)    ###### Outer (and only) loop.
  112.   {
  113.     k = k+1;
  114.  
  115.     fmax = f[1];
  116.     if (fmax > fmax_old)
  117.     {
  118.       if (!isempty(savit))
  119.       {
  120.         x = reshape (V[;1], x.nr, x.nc); 
  121.         write("savit", x, fmax, nf);
  122.       }
  123.       if (trace)
  124.       {
  125.         printf("Iter. %2.0f"', k);
  126.         printf("  how = %s  ", how);
  127.         printf("nf = %3.0f,  f = %9.4e  (%2.1f)\n", nf, fmax, ...
  128.                100*(fmax-fmax_old)/(abs(fmax_old)+eps));
  129.       }
  130.     }
  131.     fmax_old = fmax;
  132.  
  133.     ### Three stopping tests from MDSMAX.M
  134.  
  135.     # Stopping Test 1 - f reached target value?
  136.     if (fmax >= stopit[3])
  137.     {
  138.       msg = "Exceeded target...quitting\n";
  139.       break  # Quit.
  140.     }
  141.  
  142.     # Stopping Test 2 - too many f-evals?
  143.     if (nf >= stopit[3])
  144.     {
  145.       msg = "Max no. of function evaluations exceeded...quitting\n";
  146.       break  # Quit.
  147.     }
  148.  
  149.     # Stopping Test 3 - converged?   This is test (4.3) in [1].
  150.     v1 = V[;1];
  151.     size_simplex = norm(V[;2:n+1]-v1[;ones(1,n)],"1") / max([1, norm(v1,"1")]);
  152.     if (size_simplex <= tol)
  153.     {
  154.       sprintf(msg, "Simplex size %9.4e <= %9.4e...quitting\n", ...
  155.               size_simplex, tol);
  156.       break  # Quit.
  157.     }
  158.  
  159.     #  One step of the Nelder-Mead simplex algorithm
  160.     #  NJH: Altered function calls and changed CNT to NF.
  161.     #       Changed each `fr < f(1)' type test to `>' for maximization
  162.     #       and re-ordered function values after sort.
  163.  
  164.     vbar = (sum(V[;1:n]')/n)';    # Mean value
  165.     vr = (1 + alpha)*vbar - alpha*V[;n+1]; 
  166.     x = reshape (vr, x.nr, x.nc);
  167.     fr = fun (x);
  168.     nf = nf + 1;
  169.     vk = vr;  fk = fr; how = "reflect, ";
  170.     if (fr > f[n])
  171.     {
  172.       if (fr > f[1])
  173.       {
  174.         ve = gamma*vr + (1-gamma)*vbar; 
  175.         x = reshape (ve, x.nr, x.nc);
  176.         fe = fun (x);
  177.         nf = nf + 1;
  178.         if (fe > f[1])
  179.         {
  180.           vk = ve; fk = fe;
  181.           how = "expand,  ";
  182.         }
  183.       }
  184.     else
  185.       vt = V[;n+1];
  186.       ft = f[n+1];
  187.       if (fr > ft)
  188.       {
  189.         vt = vr;
  190.         ft = fr;
  191.       }
  192.       vc = beta*vt + (1-beta)*vbar; 
  193.       x = reshape (vc, x.nr, x.nc);
  194.       fc = fun (x);
  195.       nf = nf + 1;
  196.       if (fc > f[n])
  197.       {
  198.         vk = vc; fk = fc;
  199.         how = "contract,";
  200.       else
  201.         for (j in 2:n)
  202.         {
  203.           V[;j] = (V[;1] + V[;j])/2;
  204.           x = reshape (V[;j], x.nr, x.nc);
  205.           f[j] = fun (x);
  206.         }
  207.         nf = nf + n-1;
  208.         vk = (V[;1] + V[;n+1])/2; 
  209.         x = reshape (vk, x.nr, x.nc);
  210.         fk = fun (x);
  211.         nf = nf + 1;
  212.         how = "shrink,  ";
  213.       }
  214.     }
  215.     V[;n+1] = vk;
  216.     f[n+1] = fk;
  217.     temp = sort(f).val;
  218.     j = sort (f).ind;
  219.     j = j[n+1:1:-1];
  220.     f = f[j]; 
  221.     V = V[;j];
  222.  
  223.   }    ###### End of outer (and only) loop.
  224.  
  225.   # Finished.
  226.   if (trace) { printf(msg); }
  227.   x = reshape (V[;1], x.nr, x.nc);
  228.  
  229.   return <<x = x; fmax = fmax; nf = nf>>;
  230. };
  231.